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1  Introduction 


A  realistic  mathematical  model  for  a  watershed-reservoir-dam  system  should 
involve  some  inherently  random  components  leading  to  a  probabilistic  represen¬ 
tation.  A  closed-form  solution  for  the  response  of  the  system  to  a  given  water- 
input  event  proves  difficult  to  obtain,  if  not  impossible,  in  view  of  the  complexity 
and  ultimate  nonlinearity  of  the  governing  formulations.  Consequently,  to  arrive 
at  useful  conclusions  in  the  study  of  a  given  system,  the  analyst  must  conduct  a 
large  nmnber  of  numerical  simulations  of  possible  realizations.  For  this  proce¬ 
dure  to  be  practical  in  exercises  of  risk  assessment  of  spillways  and  dam 
nonoverflow  monoliths  imder  flood  hazard,  the  associated  mathematical  model 
must  be  effective  and  economical.  This  report  introduces  a  new  engineering 
model  for  flood-risk  analysis  of  watershed-reservoir-dam  systems  designed  for 
application  in  multiple  mathematical  simulations  while  capturing  the  essential 
characteristics  of  the  physical  system. 

Hyetographs  for  a  given  basin  are  built  by  inserting  rational  elements  into  the 
cxirrent  state  of  practice  (U.S.  Bureau  of  Reclamation  (USER)  1976, 1977).  A 
simplified  representation  of  the  watershed  is  introduced  at  the  level  of  the  fun¬ 
damental  unitgraph.  All  subsequent  compositions  are  mathematically  rigorous, 
leading  to  a  convolution  integral  for  a  rational  watershed-output  hydrograph. 

The  reservoir  is  represented  by  a  nonlinear  ordinary  differential  equation  formu¬ 
lated  on  the  basis  of  the  principle  of  continuity  and  on  the  assumption  of  a  con¬ 
vex  reservoir.  The  spillway  discharge  is  modeled  using  the  von  Mises 
semiempirical  expression  for  a  wide-weir  flow  (Street,  Walters,  and  Vennard 
1996).  Deterministic  interpretations  of  the  model  provide  insight  into  physical 
behavior  through  parametric  studies  on  the  occurrence  of  hydrograph  peaks,  re¬ 
sponse  spectra,  and  residual  reservoir  pools.  Stochastic  interpretations  of  the 
model  provide  insight  into  the  resulting  response  random  processes  and  the  asso¬ 
ciated  hazard  curves  necessary  for  subsequent  evaluations  of  the  probabilities  of 
overtopping,  overstressing,  overtumiug,  and  sliding  failures  (Ellingwood  1995, 
de  Bejar  1999). 

The  relevant  features  of  the  USER  recommended  practical  procedme  are  re¬ 
tained  in  the  model.  Soil  Conservation  Service  (Miller  and  Clark  1960)  maps  of 
probable  maximum  precipitation  (PMP)  during  extreme  storms  are  adopted  as 
mean  values  of  extreme-value  distributions.  USER  empirical  factors  for  the  de¬ 
termination  of  excess  rain  are  directly  implemented  as  are  the  Soil  Conservation 
Service  (SCS)  charts  for  cover-and-land-use  complex  coefficients. 
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To  focus  on  tile  effect  of  a  few  essential  factors  on  the  system  response  and 
to  keep  the  formulation  sufficiently  simple  to  promote  physical  insight,  only  the 
basic  random  variables  are  included  in  the  watershed  stochastic  model:  (a)  storm 
magnitude  and  (b)  watershed  characteristic  centroidal  lag  time.  The  study  of  the 
random  process  representing  the  variation  in  reservoir-pool  in  response  to  the 
water-input  event  is  conducted  with  and  without  the  presence  of  noise  in  the 
reservoir,  for  mutual  comparison  among  the  resulting  hazard  curves. 
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The  physical  system  to  be  modeled  consists  of  three  major  components  (Fig¬ 
ure  1):  the  watershed  or  drainage  basin,  the  reservoir  or  excess-rain  storage,  and 
the  spiUway-dam  stracture,  generally  equipped  with  a  gate  system  for  flood  con¬ 
trol  and  evacuation. 

The  U.S.  Weather  Bureau  in  collaboration  with  the  U.S.  Army  Corps  of  En¬ 
gineers  have  developed  empirical  charts  to  estimate  the  6-hr,  lO-square-mile  PMP 
as  the  result  of  a  uniform  storm  on  such  ‘a  point’  within  the  basin  area  (USER 
1976,  1977).  The  proper  chart  to  be  applied  depends  on  the  specific  geographical 
location  of  the  project  site.  The  charts  are  built  for  U.S.  zones  either  east  or  west 
of  the  105°  meridian.  To  be  specific,  the  model  developed  here  apphes  to  U.S. 
watersheds  east  of  the  105°  meridian,  but  a  parallel  development  may  just  as  eas¬ 
ily  be  formulated  for  U.S.  western  watersheds. 

Subsequently,  the  point-storm  PMP  is  scaled  up  on  the  basis  of  empirical 
charts  for  the  size  of  the  specific  drainage  area  under  consideration  and  for  several 
values  of  storm  duration  to  estimate  upper  bounds  of  cumulative  total  rain  falhng 
uniformly  over  the  watershed.  These  estimates  of  precipitation  are  taken  as  the 
known  PMP  distribution  in  time  in  the  deterministic  formulations,  or  as  mean  val¬ 
ues  of  the  extreme-value  distribution  of  the  largest  values,  type  I  (Gumbel  distri¬ 
bution)  in  the  probabilistic  formulations. 

The  watershed  represents  the  first  filter  in  the  system.  A  substantial  portion  of 
the  falling  rain  is  lost  as  the  result  of  a  variety  of  factors.  Among  the  main  factors 
contributing  to  rain  loss  are  (a)  evaporation  and  transpiration,  (b)  retention  by 
vegetation  and  by  topographic  details  of  the  terrain  (including  minor  ponds),  and 
(c)  surficial  infiltration  and  deep  percolation,  depending  on  the  type  of  soil  cover 
and  geological  characteristics  of  the  region.  The  difference  between  total  rain  and 
losses  other  than  those  from  rapid-drainage  flowing-through-cover  soil  water  is 
the  excess  rain  defined  here  as  runoff. 

Runoff  is  evacuated  relatively  rapidly  from  the  watershed  via  open  chaimels 
and  rivers  that  lead  to  the  reservoir  entrance  (point  B  in  the  schematic  representa¬ 
tion  in  Figure  1),  where  the  runoff  flow  Q(t)  can  be  measured  (Figure  2). 

The  storm  total  rain  is  translated  into  incremental  runoff  using  semi-empirical 
transformations  (Miller  and  Clark  1960)  that  include  consideration  of  a  complex 
index  to  classify  the  watershed  according  to  the  cover  soil  and  the  land  use.  This 
piecewise-constant  effective  water-input  history  is  apphed  to  the  watershed  hourly 
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Figure  1 .  Schematic  plan  view  of  a  watershed-reservoir-dam  system 


Figure  2.  Longitudinal  section  of  analytical  model  of  a  reservoir 
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during  the  first  6  hr  of  storm  and,  thereafter,  in  incremental  time  intervals  of  6, 12, 
and  24  hr,  respectively.  It  should  be  noted,  however,  that  the  order  of  hourly  pre¬ 
cipitation  during  the  fct  6  hr  of  storm  cannot  be  predicted.  Therefore,  in  this 
model  the  specific  order  of  the  rain  steps  for  a  given  storm  is  subjected  to  aleatory 
permutation,  whether  the  representation  is  deterministic  or  not. 

Figure  2  shows  the  longitudinal  section  of  the  analytical  model  for  the  reser¬ 
voir  component  in  the  system.  The  input  flow  Q(t)  represents  the  response  of  flie 
watershed  to  the  water-input  event.  The  spillway  outflow  q{h(t),  c)  at  the  opposite 
end  of  the  reservoir  (Figure  3)  depends  on  the  discharge  coefficient  c.  This,  in 
turn,  is  a  function  of  the  weir  elevation  z  and  of  the  reservoir  elevation  h(t)  itself, 
giving  rise  to  a  highly  nonMnear  governing  differential  equation  for  the  response. 
The  reservoir  water  level  over  the  spillway  crest  h(t)  is  directly  related  to  the 
storm-water  storage  and  represents  the  hazard  on  the  dam  structure  whose  safety 
is  to  be  evaluated  later  in  terms  of  the  potential  modes  of  failure  of  overtopping, 
overstressing,  overturning,  and  sUding  instabihty.  For  simplicity,  the  pool  level  at 
the  beginning  of  the  storm  is  considered  to  be  that  of  the  spillway  crest,  but  any 
other  convenient  datum  may  be  defined  by  the  analyst  at  wiU. 


Figure  3.  Transverse  section  of  analytical  model  on  nonoverflow  dam  monolith 
(Spillway  crest  is  schematically  represented  by  the  dotted  line) 
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Notice  that  the  model  differs  from  the  real  reservoir  by  the  fact  that  the  bot¬ 
tom  of  the  reservoir  may  be  randomly  covered  with  sediments  over  time  and  also 
that  our  topographical  surveys  are  imperfect,  particularly  in  the  vicinity  of  the 
boundary  of  the  reservoir.  The  corresponding  effects  on  the  governing  differential 
equation  for  the  reservoir  response  are  optionally  included  in  the  probabilistic 
version  of  the  model  by  adding  a  Gaussian  noise  component  to  the  input  random 
process. 


Convex  Unitgraph 

The  convex  model  for  the  response  of  a  watershed  to  a  water-input  event  is 
based  on  the  principle  of  continuity  and  on  a  postulated  linear  relationship  be¬ 
tween  the  response  flow  Q(t)  and  the  watershed  storage  S{t).  When  the  storm  in¬ 
flow  is  constant,  the  governing  differential  equation  may  be  expressed  (Dingman 
1994)  as: 


dS 
—  + 
dt 


5  =  6)0 


(1) 


where  T*  =  centroidal  lag  between  the  inflow  and  outflow  hydrographs 

This  is  taken  as  a  constant  characteristic  of  the  watershed.  The  solution  to  this 
equation  under  zero  initial  conditions  is: 


where  k*  -  l/T* 

This  result  indicates  that  the  output  response  approaches  asymptotically  the  inflow 
cOo  as  far  as  the  termination  of  the  water-input  event  toy  At  this  time,  the  recession 
limb  of  the  response  starts  an  exponential  decay  toward  the  zero-flow  value.  This 
means  that  the  time  to  concentration  of  the  model  tends  to  infinity,  similarly  to  the 
time  to  concentration  of  a  real  watershed.  The  recession  limb  of  the  response  is 
given  by: 

6(0  =  qpk-  ^  t>to)  (3) 


where  qpk  =  Q(ta,)  is  the  maximum  value  of  the  outflow 

Notice  that  the  only  parameter  characterizing  the  response  is  T*,  which  may  be 
estimated  as: 


y*  = 


h. 

u^, 


(4) 


where 

Lr  =  length  of  the  longest  reach  in  the  drainage  basin  (distance  between 
points  A  and  B  along  the  stream  in  Figure  1) 
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Uw  =  velocity  of  propagation  of  the  flood  wave 

The  convex  unitgraph  is  defined  in  this  investigation  as  the  outflow  response 
of  the  watershed  to  a  constant-rate  water-input  event  with  a  total  volume  of  unit 
value  (usually  a  volume  of  rain  with  1  in.  of  depth  and  uniformly  distributed  over 
the  whole  drainage  area).  The  duration  of  this  constant-rate  water-input  event  is 
the  unit  period,  which  in  this  model  is  taken  as  1  hr. 


Watershed  Response  Hydrograph 

The  watershed  response  hydrograph  is  obtained  by  superposition  in  time  of 
the  scaled  hydrographs  corresponding  to  the  actual  incremental  effective  rain  vol¬ 
umes  (usually  expressed  in  inches  of  rain  uniformly  distributed  over  the  whole 
watershed  area).  Analytically,  this  procedure  can  be  generalized  by  considering 
the  response  to  a  unit-volume  input  rain  concentrated  at  the  origin  of  time.  In 
other  words,  the  watershed  output  response  is,  in  this  case,  the  unit-impulse  re¬ 
sponse  function  u(t),  as  the  input  rain  is  a  Dirac  delta  function  at  the  origin  of 
time,  i.e.,  a  zero-duration  unit-volume  rain  at  r  =  0.  In  this  case,  the  governing 
differential  equation  for  the  watershed  becomes; 

T*.^+u{t)  =  8{t)  (5) 

dt 

with  initial  condition  u(0)  =  0.  The  general  solution  of  this  equation  is: 

u{t)  =  Ae-^*-^  (6) 

The  constant  of  integration  A  may  be  obtained  by  integrating  Equation  5  over 
the  infinitesimal  time  interval  (-e,+e)  and  taking  the  limit  as  e->0  to  get  the  unit- 
impulse  response  function  as: 

=  f>0  (7) 

The  watershed  response  hydrograph  to  an  inflow  with  the  rate  of  rain  r(t)  may 
be  obtained  by  superposition  in  the  time  domain  and  is  given  by  the  convolution 
integral: 


Q[t)=  j  r{t-T)  u{'v)  dT:  t>0 

0+ 


(8) 


In  fact,  the  response  to  r(t)  =  (0^=  constant  may  also  be  obtained  as  the  expres¬ 
sions  in  Equations  2  and  3  by  direct  evaluation  of  Equation  8. 


Figure  4  shows  a  family  of  these  responses  for  a  storm  with  magnitude  tab 
=  100  Ksec-ft  acting  during  4  hr.  The  parameter  T*  for  the  family  of  curves  var¬ 
ies  from  a  value  T*  =  25  hr  (a  slow-evacuation  watershed)  to  a  value  T*  =  2.5  hr 
(a  rapid-evacuation  watershed).  The  ‘slow’  watersheds  have  smaller  peak  values 
of  outflow  and  larger  residual  flows  at  the  end  of  the  period  under  consideration 
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Figure  4.  Family  of  convex-watershed  responses  to  a  uniform  runoff  flow  Wq 
acting  during  a  finite  time  interval  (4  hr).  Curves  in  the  set  are 
characterized  by  the  centroidal  lag  time  T* 

(10  hr).  On  the  other  hand,  the  ‘rapid’  watersheds  often  get  almost  to  the  asymp¬ 
totic  value  tBb  within  the  storm  duration  and  rapidly  decay  toward  zero  residual 
flow  upon  storm  termination. 

The  procedure  described  up  to  this  point  to  construct  the  inflow  design  flood 
hydrograph  into  the  reservoir  (i.e.,  the  watershed  outflow  hydrograph)  is  shown  in 
the  flow  diagram  in  Figure  5.  This  schematic  stracture  is  followed  in  both  the 
deterministic  and  the  probabilistic  versions  of  the  model  being  constructed. 

Again,  the  probabihstic  model  concentrates  on  the  effects  of  a  random  water-input 
event  magnitude  and  of  a  random  watershed  centroidal  lag  T*. 
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Figure  5.  Procedure  to  construct  an  inflow-flood  hydrograph  for  a  typical 
reservoir 
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3  Deterministic  Studies 


An  examination  of  the  physical  quantities  entering  the  detenninistic  formula¬ 
tion  for  the  response  of  the  watershed-reservoir-dam  system  to  a  water-input  event 
and  its  sensitivity  to  variations  in  those  quantities  proves  to  be  insightful  in  re¬ 
vealing  the  fundamental  nature  of  the  relations  involved. 


Watershed  Routing 

In  this  model,  the  response  of  the  watershed-reservoir-dam  system  is  followed 
for  a  time  interval  of  60  hr  after  the  begiiming  of  the  storm,  which  is  assumed  to 
last  48  hr.  The  subscript  I  =  1,2,. ..,6  refers  to  the  first  six  l-hr  intervals;  I  =  7,8,9 
refers  to  the  time  intervals  following  the  first  6  hr  of  storm  duration  with  time 
lengths  of  6,  12,  and  24  hr,  respectively;  and  i  =  10  refers  to  the  last  12  hr  of  rain- 
free  history  of  response. 

The  procedure  to  determine  the  incremental  runoff  during  each  of  these  inter¬ 
vals  follows  the  standard  and  well  accepted  recommendations  in  practice  (Miller 
and  Clark  1960;  USSR  1976, 1977).  According  to  the  geographical  location  and 
extension  of  the  watershed,  the  6-hr  ‘point’  PMP  is  distributed  in  time  and  ex¬ 
pressed  as  the  accumulative  rain  fall  in  the  actual-size  watershed  at  the  end  of 
each  of  the  intervals  described  above  (a  ,  i  =  6,..,9).  This  distribution  may  be  ex¬ 
pressed  as: 

n+5=aiP  i  =  l,...,4  (9) 


where 

p  =  6-hr  ‘point’  PMP 
a,  =  empirical  coefficient 

Likewise,  the  accumulative  rainfall  at  the  end  of  each  of  the  first  six  l-hr  intervals 
may  be  expressed  as: 

n=bir^  i=l,...,6  (10) 

where 

rg  =  the  accumulative  total  rain  at  the  end  of  6  hr  of  storm,  as  provided  by 
Equation  9 


10 


Chapter  3  Deterministic  Studies 


bj  =  empirical  coefficient 

Therefore,  the  incremental  total  rain  corresponding  to  each  of  the  first  1-hr  inter¬ 
vals  is  given  by: 

i  =  2,...,6 

Again,  the  order  of  these  first  six  1-hr  incremental  contributions  to  the  total 
precipitation  cannot  be  predicted,  and  they  are  given  an  aleatory  permutation,  after 
which  the  accumulative  total  rain  is  recalculated  according  to: 


1  =  2,.. .,9 


(12) 


Next,  the  accumulative  and  incremental  effective  precipitations  (direct  runoff) 
are  calculated.  These  calculations  require  the  estimation  of  the  local  hydrologic 
soil-cover  complex  number  (S),  according  to  the  soil-series  classification  and  the 
combined  land  use  at  the  site  (USER  1976, 1977),  modified  according  to  the  an¬ 
tecedent  conditions.  The  U.S.  Soil  Conservation  Service  (Miller  and  Clark  1960) 
recommends  the  use  of  the  following  fit  to  estimate  runoff  (based  on  numerous 
statistical  studies  and  assuming  the  initial  abstraction  as  /«  =  0.2-S): 


{n-o.2-sf 

ri-O.SS 


i  =  l,...,9 


where  p,  =  accumulative  runoff  at  the  end  of  the  i-th  time  interval,  and 

Api=Pi 

^i=Pi-Pi-l  i  =  2,...,9 


(13) 


(14) 


where  Api  =  incremental  runoff  corresponding  to  the  i-th  time  interval 

There  is  a  physical  lower  bound  for  the  hourly  loss  during  each  interval 
(Miller  and  Clarkl960).  In  this  model,  the  hourly  loss  is  not  allowed  to  be  less 
than  0.05  in.  Upon  insertion  of  this  miitimum  value,  the  incremental  runoff  is 
recalculated  for  each  interval  as: 


Api=Ari-ALi  i  =  l,...,9 


(15) 


and 


*  * 

ft  =  Aft 

Pi=Ap*  +  Pi_i  I  =2,..., 9 

where  ALj  is  the  total  loss  corresponding  to  the  i-th  interval,  and  symbols  with  an 
attached  asterisk  correspond  to  the  recalculated  quantities. 
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The  watershed  routing  is  completed  using  the  convolution  integral  of  Equa¬ 
tion  8  reformulated  for  computational  analysis.  When  the  input  effective  precipi¬ 
tation  is  represented  by  a  piecewise  constant  function,  Equation  8  gives 

2(0  =  '^<lpk,i  •  t  >  to, 

i=\ 


2(0  =  '^{Recess). 

i=l 


(17) 


where 

(Recess)i  =  contribution  of  the  recession  limb  of  the  response  to  the  i-th  in¬ 
crement  of  direct  runoff  when  acting  alone  dining  the  time  inter¬ 
val  ending  at  t, 

qpki  =  corresponding  maximum  response,  as  in  Equation  3,  and 


2(0 


*j<t<ta  j^O 


(18) 


when  the  response  is  evaluated  within  the  (/+l)-th  interval  and  to  =  0 


Reservoir  Routing 


The  outflow  from  the  routing  of  rain  fall  through  the  watershed  represents  the 
inflow  for  the  reservoir  component.  Based  on  the  principle  of  continuity,  the  gov¬ 
erning  differential  equation  for  the  reservoir  routing  in  terms  of  the  pool  at  the 
upstream  face  of  the  dam,  h(t),  is  given  by  (Jiang  1998): 


dh  _  Q{t)-q[h{t),c 
dt  ~  G\h{t)\ 


(19) 


where 

Q( t)  =  inflow  from  the  watershed  given  by  Equations  1 7  and  1 8 

q(h( t),c)  =  outflow  through  the  spillway  weir 

c  =  corresponding  discharge  coefficient 

G(h)  =  gradient  of  variation  of  the  reservoir  storage  as  a  function  of  the 
reservoir  pool 

Notice  that  other  outlet  works  may  be  included  in  q(h(t),c),  if  present.  In  this 
model,  for  simplicity,  only  a  rectangular  spillway  weir  is  considered.  The  outflow 
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through  a  rectangular  weir  may  be  estimated  by  (Street,  Watters,  and  Vennard 
1996): 

q[h{t),c'\  =  ^-bc-yf2g (20) 


where 

b  =  weir  length 
c  =  discharge  coefficient 
g  =  acceleration  of  gravity 

Von  Mises  has  developed  a  simple  semiempirical  expression  for  the  coefficient  of 
discharge  as  (Olson  1961): 

c  =  0.61 1  +  0.075-^  (21) 

z 

where  z  is  the  height  of  the  spillway  crest  over  the  reservoir  bottom  (Figure  3). 

The  reservoir  storage  may  be  regressed  empirically  as  a  quadratic  function  of 
the  reservoir  pool  (Jiang  1998),  leading  to  a  linear  fit  for  the  corresponding  gradi¬ 
ent,  so  that  it  may  be  expressed  as: 

G{h)  =  a  +  ph  (22) 

where  are  empirical  coefficients  for  the  reservoir  being  considered. 

Inserting  Equations  17, 18,  20,  21,  and  22  into  Equation  19  leads  to  a  nonlinear 
ordinary  differential  equation  for  h(t)  that  can  be  solved  numerically  (The  Math- 
Works  2000). 


Computational  Analysis 

The  watershed  centroidal  lag  is  estimated  in  the  computational  implementa¬ 
tion  of  the  model  using  the  empirical  equation  inferred  by  the  CaUfomia  High¬ 
ways  and  Public  Works  (Miller  and  Clark  1960): 


T*  = 


/  T  nO.SBS 

11.9L^ 


H 


(23) 


where 

L  =  length  of  the  longest  watercourse  in  miles  (i.e.,  the  distance  between 
points  A  and  B  along  the  stream  in  Figure  1) 

H  =  elevation  difference  (between  the  same  points  A  and  B  in  Figure  1) 
in  feet 
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T*  =  obtained  in  hours 

Figure  6  shows  a  Graphical  User  Interface  (GUI)  implemented  in  Matlab  6.0 
(The  MathWorks  2000)  with  the  data  for  a  sample  watershed-reservoir-dam  sys¬ 
tem  in  zone  7  east  of  the  105°  meridian.  The  program  builds  the  hydrographs  for 
the  inflow  design  flood  into  the  reservoir  and  for  the  pool,  for  the  design  storm 
(PMP),  and  provides  the  analyst  with  instant  feedback.  The  command  Plot  Hy¬ 
drographs  on  the  lower  right  comer  of  the  GUI  invokes  the  callback  routines  that 
execute  the  convex  routing  algorithm  with  immediate  results  useful  for  either  the 
design  of  new  hydrauUc  structures  or  the  evaluation  of  existing  facilities  against 
overtopping. 


Watershed  Outflow  Hydrograph 
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Reservoir  Pool  Hydrograph 
25  r 


0*- 

0 


20  40 


60 


time  [hr] 


regress  param 

|_  pjipt  t-ydrographs  ] 


iSSTsTi? 


I  beta  -  ’  1E*2 


Figure  6.  GUI  designed  for  ready  estimation  of  deterministic  design 
hydrographs 


Surface  Hydrographs 

The  watershed  outflow  hydrograph  and  the  associated  reservoir  pool  hydro¬ 
graph  are  then  generated  for  a  drainage  basin  with  a  continuously  varying 
centroidal  lag.  Figures  7  and  8  show  the  resulting  siuface  hydrographs  for  a 
rather  severe  design  storm  with  PMP  =  30  in.  Notice  the  pronounced  peak  of  the 
hydrographs  (i.e.,  of  the  trace  of  the  corresponding  surface  on  the  plane  T*  = 
constant)  for  watersheds  with  rapid  evacuation  characteristics  (relatively  short  T*) 
and  the  small  residual  quantity  at  the  end  of  the  response  history.  Watersheds 
with  slow  evacuation  characteristics  (relatively  long  T*)  tend  to  remain  flat  in 
time,  and  therefore  show  larger  residual  quantities  at  the  end  of  the  response 


14 


Chapter  3  Deterministic  Studies 


history.  Figures  9,  10,  11,  and  12  show  these  hydrographs  for  selected  PMP  spe¬ 
cific  values  of  30  and  20  in.,  respectively.  The  trends  in  the  family  of  curves  re¬ 
main  the  same  with  the  smaller  magnitudes  of  response  associated  with  the 
smaller  design  storms. 

Similarly,  hydrographs  for  watershed  outflow  and  for  the  corresponding  res¬ 
ervoir  pool  are  built  for  a  fixed  watershed,  with  a  varying  magnitude  of  the  water- 
input  event  (PMP  =  15  to  35  in.,  @5  in.).  Figures  13  and  14,  and  15  and  16  show 
two  sets  of  hydrograph  families  for  a  rapid-evacuation  watershed  (T*  =  5  hr)  and 
for  a  slow-evacuation  watershed  (T*  =  15  hr).  Both  sets  tend  to  support  the  basic 
assumption  under  the  standard  unitgraph  superposition  principle  applied  in  prac¬ 
tice  that  the  response  hydrographs  essentially  retain  their  shapes  with  the  ordinates 
amplified  in  proportion  to  the  magnitude  of  the  water-input  event  (PMP). 
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RESERVOIR  POOL  HYDROGRAPH  (T*=5hr) 


Figure  14.  Family  of  reservoir-pool  hydrographs  built  for  a  deterministic  v/aXer- 
shed  with  a  relatively  short  centroidal  lag  time  (7*=  5  hr)  for  several 
values  of  the  water-input  event  PMP 


Figure  15.  Family  of  reservoir  inflow-flood  hydrographs  built  for  a  deterministic 
watershed  with  a  relatively  long  centroidal  lag  time  (7*  =  15  hr)  for 
several  values  of  the  water-input  event  PMP 
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RESERVOIR  POOL  HYDROGRAPH  (r=15hr) 


Figure  16.  Family  of  reservoir-pool  hydrographs  built  for  a  deterministic  water¬ 
shed  with  a  relatively  long  centroidal  lag  time  (T*  =  15  hr)  for  several 
values  of  the  water-input  event  PMP 

Response  Spectra 

The  peak  values  of  response  in  terms  of  the  watershed  outflow  and  the  reser¬ 
voir  pool  may  be  recorded  for  a  continuous  variation  of  the  centroidal  lag  and  with 
the  magnitude  of  the  water-input  event  (PMP)  as  a  parameter.  The  resulting  fam¬ 
ily  of  curves  constitute  the  response  spectra  for  the  response  quantity  under  con¬ 
sideration.  Figures  17  and  18  show  the  response  spectra  for  the  watershed 
outflow  Q(t)  and  for  the  reservoir  pool  h(t). 

The  set  of  response  spectra  for  the  watershed  outflow  is  very  sensitive  to  the 
particular  order  of  the  resulting  aleatoiy  permutation  of  the  hourly  storm  precipi¬ 
tation  values  during  the  first  6  hr  of  the  water-input  event.  More  unfavorable 
system  responses  are  obtained  when  the  largest  incremental  rain  value  occurs 
during  the  sixth  storm  hour  than  the  responses  for  those  storms  in  which  the  larg¬ 
est  incremental  rain  value  occurs  during  the  first  hour.  This  sensitivity  becomes 
apparent  for  the  rapid-evacuation  watersheds  (those  with  relatively  short  T*),  as 
evidenced  by  the  nonsmooth  gradient  of  Q,max  with  respect  to  the  storm  magni¬ 
tude  (PMP)  in  Figure  17. 

By  contrast,  the  set  of  spectra  for  the  reservoir  pool  exhibit  relatively  smooth 
transitions  when  the  storm  magnitude  (PMP)  is  varied.  There  are  two  fundamen¬ 
tal  reasons  for  this  behavior:  (a)  the  peak  response  of  h(t)  for  a  given  hydrograph 
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systematically  shows  a  delay  with  respect  to  the  peak  of  Q(t)  (Figure  6),  occurring 
at  some  time-distance  of  the  particular  hourly  incremental  rain  lumps  during  the 
first  6  hr  of  storm,  and  (b)  the  peak  response  of  h(t)  occurs  after  a  second  filter  in 
the  system  has  acted  on  the  water  input  to  the  watershed.  The  reservoir  pool 
spectra  in  Figure  18  can  be  closely  approximated  by  the  following  useful  mathe¬ 
matical  expression: 

hm^=h2{T*)  +  ^  ^  ^  ^-^  {PMP-IS)  PMP>15  (24) 


where 

/H  (r  *)  =  47.15- 8.416 -^(7*) 

hi{T*)  =  16.82  -  0.5192  •  T  *40.007423 (T  *f 
and  the  PMP  is  given  in  inches  of  rain. 


Residual  Response 

Similarly,  the  residual  values  of  response  in  terms  of  the  watershed  outflow 
and  the  reservoir  pool  at  the  end  of  the  observed  history  may  be  recorded  and  rep¬ 
resented  graphically,  as  in  Figures  19  and  20,  respectively.  Now  the  recorded 
values  take  place  at  a  far  time-distance  of  the  local  rain  distribution  dming  the  first 
6  hr  of  storm,  and  the  resulting  curves  exhibit  smooth  gradients  with  the  magni¬ 
tude  of  the  water-input  event  (PMP).  The  watershed  residual  outflow  in  Figure  19 
can  be  closely  approximated  by  the  following  mathematical  expression: 

Qr  =Qr2{T*)+  B  . ^pMP - 1 5)  PMP >  15  (25) 


where 

0;.!  (r  *)  =  -0.6571  +  0. 1 253  •  r  *  -K).01 648  {T*f  -  0.0003444  •  (T 

Qrl  {T  *)  =  -0.06143  -  0.01471  •  T  *  40.008029(7’  *f  -  0.0001489  (T  *f 
and  the  PMP  is  given  in  inches  of  rain. 
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WATERSHED  OUTFLOW  RESIDUAL  @  1=60  hr 


Figure  19.  Family  of  residual  watershed  outflows,  at  the  end  of  the  history  being 
considered  (t  =  60  hr),  for  several  values  of  the  water-input  event  PMP 


Figure  20.  Family  of  residual  reservoir  pools,  at  the  end  of  the  history  being 

considered  (t  =  60  hr),  for  several  values  of  the  water-input  event  PMP 
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Likewise,  the  reservoir  residual  pool  in  Figure  20  can  be  closely  approximated 
by  the  following  mathematical  expression: 

/i^=/V2(r*)+— ^^-(PMP-IS)  PMP>\5  (26) 

where 

K\  {T  *)  =  3.481  +  0.3569  •  T  *  +0.001561  (r  *f  -  0.0001452  •  {T  *f 
and 

h^2  [T  *)  =  1.921  +  0.2061  •  T  * -0.001773(7’  *f  -  {1.2,11  E  - 5) •  (r  *f 
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4  Probabilistic  Studies 


Instantaneous  Response  Distributions 

Next,  two  fundamental  variables  into  the  routing  of  the  water-input  event 
through  the  watershed-reservoir-dam  system,  namely  the  storm  magnitude  and  the 
watershed  centroidal  lag,  are  assumed  to  be  random.  The  6-hr  ‘point’  storm  mag¬ 
nitude  is  modeled  with  the  extreme-value  distribution  type  I,  of  the  largest  values 
(Gumbel  distribution),  with  mean  value  given  by  the  deterministic  assessment  of 
the  6-hr  ‘point’  PMP  (MiUer  and  Clark  1960;  USER  1976, 1977)  and  with  a  coef¬ 
ficient  of  variation  estimated  as  0.1.  The  watershed  centroidal  lag  is  assumed  as 
log-normally  distributed,  with  mean  value  given  by  the  deterministic  assessment 
and  with  a  coefficient  of  variation  estimated  as  0.3.  Our  objective  in  this  investi¬ 
gation  is  to  identify  the  instantaneous  distribution  of  the  response  hydrograph  for 
the  reservoir  pool,  h(t),  since  this  represents  the  effective  hazard  on  the  dam 
(de  Bejar  1999).  Sampling  is  conducted  on  the  basis  of  Montecarlo  simulations 
(Benjamin  and  Cornell  1970). 

The  best  distribution  is  judged  to  be  the  log-normal  distribution,  fitting  the 
data  well,  while  being  sufficiently  simple  for  engineering  applications.  Figure  21 
shows  the  quantile-quantile  plot  for  the  random  variable  natural  logarithm  of  the 
reservoir  pool  at  time  =  24  hr,  i.e.,  ln[h(24)],  resulting  from  10,000  simulated  re¬ 
alizations  of  the  routing  of  a  random  storm  through  a  random  watershed,  with  aU 
other  parameters  in  the  deterministic  system  described  in  Figure  6.  The  approxi¬ 
mate  alignment  with  a  straight  line  is  indicative  of  the  goodness-of-fit  for  engi¬ 
neering  applications  of  the  underlying  log-normal  assumption  for  the  distribution 
of  h(24).  Figure  22  shows  the  nonparametric  probability  density  function  corre¬ 
sponding  to  the  10,000  simulations.  Likewise,  Figures  23  and  24  show  the  corre¬ 
sponding  results  of  the  10,000  simulations  for  h(12).  Again,  the  log-normal 
distribution  is  considered  satisfactory  for  most  of  the  data,  except  possibly  at  the 
tails  of  the  distribution,  where  the  major  deviations  fi'om  the  model  occur,  as 
expected. 

Finally,  the  peak  values  of  h(t)  resulting  from  the  10,000  simulations  are  also 
processed,  and  again  the  log-normal  distribution  is  selected  as  the  most  conven¬ 
ient  for  engineering  applications.  Figure  25  shows  a  box  plot  for  the  random  vari¬ 
able  ln[h,max]  on  die  basis  of  1,000  simulations.  Only  six  deviant  points  are 
identified  at  the  upper  end  of  the  presumed  log-normal  distribution,  and  only  one 
deviant  point  at  the  lower  end  of  tibe  distribution.  In  fact.  Figure  26  shows  the 
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Ln(24-hr  pool  level  h) 


Figure  21 .  Quantile-quantile  plot  for  loge  [h(24)]  based  on  1 0,000  simulations 


Figure  22.  Nonparametric  probability  density  for  loge  [h(24)]  based  on 
1 0,000  simulations 
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Ln(12-hrpool  level  h) 
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Figure  24.  Nonparametric  probability  density  for  lope  [h(12)]  based  on 
10,000  simulations 
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Figure  25.  Box  plot  for  logs  [h.maj  based  on  1 ,000  simulations 


Normal  Distribution 


Figure  26.  Quantile-quantile  plot  for  logs  [h.maj  based  on  1 ,000 
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corresponding  quantile-quantile  plot  that  confirms  the  goodness-of-fit  of  the  log¬ 
normal  distribution  for  most  of  the  data  (again,  a  few  deviant  points  are  identified 
at  the  tails  of  the  distribution).  Figure  27  shows  the  corresponding  nonparametric 
probability  density  function  on  the  basis  of  1,000  simulations. 


Figure  27.  Nonparametric  probability  density  for  log^  [h.^ax]  based  on 
1 ,000  simulations 


Time  to  Peak  Reservoir  Pool 

The  time  (t )  to  reach  the  maximum  reservoir  pool  (h,max)  is  another  random 
variable  of  interest  in  the  routing  process.  On  the  basis  of  10,000  simulations,  the 
resultant  time  for  the  maximum  reservoir  pool  appears  as  a  cloud  of  points  with 
no  particular  pattern.  However,  upon  transformation  of  both  random  variables 
(centering  about  the  mean,  followed  by  scaling  down  by  the  corresponding  stan¬ 
dard  deviation),  the  data  are  fitted  reasonably  well  by  a  cubic  polynomial.  In 
terms  of  normalized  variables,  one  obtains: 

y  =  -0. 19  -  0.7  ■  jc  +  0.21  •  -  0.028  •  (27) 

where  y  and  x  are  the  normalized  time  and  maximum  reservoir  pool,  respectively. 
This  fit  is  shown  graphically  in  Figure  28.  Notice  that  the  sum  of  residuals 
(=  68.8)  for  this  particular  set  of  simulations  is  satisfactorily  small  for  practical 
applications.  (This  sum  of  residuals  is  not  significantly  reduced  by  increasing  the 
degree  of  the  fitting  polynomial  up  to  10  deg). 
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NORMALIZED  FIT 
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Figure  28.  Cubic  regression  to  estimate  the  time  of  occurrence  of  the  maximum 
reservoir  pool,  h.max.  based  on  10,000  simulations 

Noise  in  the  Reservoir 

Consider  that  the  bottom  of  the  reservoir  may  be  randomly  covered  with 
sediments  over  time;  that  our  topographical  surveys  are  imperfect,  particularly  in 
the  vicinity  of  the  boundary  of  the  reservoir;  or  that  the  actual  values  of  the  regres¬ 
sion  coefficients  a  and  P  for  the  gradient  of  reservoir  storage  with  changes  in  res¬ 
ervoir  pool  may  deviate  randomly  from  the  theoretical  estimations.  The  influence 
of  these  factors  on  the  reservoir  response  can  be  modeled  by  including  a  Gaussian 
noise  term  in  the  governing  differential  equation  for  the  reservoir  (Jiang  1998). 
The  result  of  such  formulation  is  equivalent  to  adding  to  the  reservoir  pool  re¬ 
sponse  obtained  before  (upon  solving  Equation  19)  the  contribution  to  the  re¬ 
sponse  from  the  Gaussian  noise  term,  h(t),  which  is  governed  by  the  following 
zero-initial  condition  Ito  stochastic  differential  equation  (Larson  and  Shubert 
1979): 


dn{t)  _  a  dWo{t) 
dt  ~  G[h{t)]  dt 
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where 


h(t)  =  total  reservoir  pool  response  (including  the  contribution  from  the 
noise) 

dWo(t)/dt  =  white  Gaussian  noise 
(T=  noise  intensity 

This  expression  may  also  be  written  as  the  following  special  case  of  the 
Langevin’s  equation: 


(29) 


where  Wq  is  the  standard  Wiener  process. 

The  solution  for  this  unpredictable  noise  is  given  by  the  following  Wiener 
integral: 


n{t)=c- 


(30) 


which  may  be  expressed  in  discrete  form  for  numerical  calculations  as  the  fol¬ 
lowing  Ito  summation: 


n{ti)  =  a 


1 

a  +  ph{tk) 


•AWo(ffc) 


(31) 


where 

AWo(4)  =  Wb(4+i) 

Wo(4)  =  zero-mean  Gaussian  process  with  variance  A4=  4+1  -  4=  At,  and 
At  =  time-step 

The  empirical  value  of  the  noise  intensity,  <J,  needs  to  be  identified  for  the 
specific  reservoir  under  study,  and  it  is  arbitrarily  taken  as  (T=  3.5E44  ft^sec'*^ 
for  numerical  calculations  in  this  investigation. 

Figure  29  shows  the  mean  function  for  the  reservoir  pool  assumed  to  be  given 
by  the  deterministic  response  hydrograph  for  the  watershed-reservoir-dam  system 
with  the  input  data  described  in  Figure  6.  The  evolution  of  the  corresponding 
response  random  process  assumed  as  log-normally  distributed  is  also  shown  in 
Figure  29.  Notice  that,  generally,  as  time  increases  along  the  recession  hmb,  there 
is  less  dispersion  in  the  prediction,  since  the  response  variances  decrease.  The 
pool  response  to  the  water-input  event  with  noise  in  the  reservoir  is  included  in 
the  figure,  for  comparison  with  the  corresponding  response  without  noise.  In  gen¬ 
eral,  the  presence  of  the  noise  increases  the  dispersion  of  the  response  random 
process  under  consideration.  Figure  30  shows  the  flood  hazard  curve  for  the  res¬ 
ervoir  pool,  with  and  without  noise  in  the  reservoir,  assuming  the  response  ran¬ 
dom  process  to  be  log-normally  distributed.  Notice  that  the  noise  causes  the 
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PROBABILITY  OF  EXCEEDENCE 


Figure  29.  Mean  function  of  the  reservoir  pool  and  schematic  representation  of 
the  random  process  h(t),  assumed  log-normal,  with  and  without  noise 
in  the  reservoir  (for  clarity,  probability  density  functions  appear 
amplified  by  a  factor  of  50) 


Figure  30.  Flood  hazard  curve  for  the  random  process  h(t),  assumed  log-normal, 
with  and  without  noise  in  the  reservoir 
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hazard  to  increase  substantially.  For  example,  given  that  the  crest  of  the  nonover¬ 
flow  section  is,  say,  30  ft  above  the  level  of  the  spillway  crest,  the  probability  of 
overtopping  turns  out  to  be  6.2  percent  for  a  noisy  reservoir,  as  compared  to 
2.4  percent  for  an  ideal  reservoir.  Figure  31  represents  a  parallel  development  for 
the  flood  hazard  curve,  assuming  that  the  reservoir  pool  response  process  is  nor¬ 
mally  distributed.  Now  the  probability  of  overtopping  is  4.5  percent  for  a  noisy 
reservoir,  as  compared  to  0.4  percent  for  an  ideal  reservoir. 
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Figure  31 .  Flood  hazard  curve  for  the  random  process  h(t),  assumed  normal, 
with  and  without  noise  in  the  reservoir 


In  general,  the  magnitude  of  the  hazard  is  increased  when  the  response  proc¬ 
esses  as  taken  as  log-normally  distributed,  as  compared  to  the  results  with  the 
normal  assumption.  The  log-normal  assumption  for  the  response  process  has  been 
proved  in  this  investigation  to  be  applicable  to  the  probabilistic  simulation  of  rout¬ 
ing  of  a  water-input  event  through  a  watershed-reservoir-dam  system.  Therefore, 
the  use  of  the  log-normal  distribution  is  recommended  for  practical  applications  to 
real  systems  in  the  field. 
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5 


Conclusions 


The  following  conclusions  are  derived  from  this  investigation: 

a.  A  rational  theoretical  model  has  been  developed  to  represent  the  routing 
of  a  water-input  event  through  a  watershed-reservoir-dam  system  and  to 
assess  its  response  in  terms  of  the  inflow  design  flood  into  the  reservoir 
and  the  resulting  reservoir  pool. 

b.  Both  deterministic  and  probabilistic  implementations  of  the  model  allow 
ready  computational  analysis,  useful  for  both  design  and  for  situational 
assessment. 

c.  System  response  random  processes  are  best  modeled  as  log-normally 
distributed. 

d  Noise  in  the  reservoir  component  is  easily  included  in  the  probabilistic 
formulation  of  the  model  and  has  an  important  influence  in  the  magnitude 
of  the  resulting  hazard.  An  accurate  assessment  of  the  noise  intensity  may 
prove  difficult  in  practice,  but  realistic  reliability  assessments  can  be  con¬ 
ducted  using  bounds  on  the  noise  intensity. 

e.  The  model  provides  direct  evaluation  of  the  probability  of  overtopping  in 
a  given  flood  scenario  and  provides  the  fundamental  hazard  curve  for 
complete  risk  analysis  of  the  dam  structure,  including  the  possible  sliding, 
overstressing,  and  overturning  modes  of  failure. 
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